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We develop the analysis of x-ray intensity correlations from dilute ensembles of 
identical particles in a number of ways. First, we show that the 3D particle struc- 
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ture can be determined if the particles can be aligned with respect to a single axis 
having a known angle with respect to the incident beam. Second, we clarify the 
phase problem in this setting and introduce a data reduction scheme that assesses 
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q ■ the integrity of the data even before the particle reconstruction is attempted. Finally, 

" c/5 . we describe an algorithm that reconstructs intensity and particle density simultane- 

^»" ously, thereby making maximal use of the available constraints. 
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1 Introduction 
> 
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As first pointed out by Kam [ 1 ], intensity fluctuations in an x-ray experiment with ensembles of non-oriented 
identical particles can provide structural information about the particles themselves. Correlations in the 
fluctuations arise from sets of photons that are scattered from the same particle and as such effectively 
provide a signal from individual particles. This technique does not place demands on the spatial coherence 
of the source — provided it spans at least one particle — but requires short pulses of high fluence and was 
largely not developed for that reason. With the availability of high intensity free-electron laser sources in 
recent years the situation has changed, and there is renewed interest in developing Kam's idea into a practical 
method. This effort has been led by Dilano and co-workers, with theoretical contributions ElO as well as 



proof-of -principle experiments |U[5). 

There is currently no known method for extracting 3D structure from intensity fluctuation data. More- 
over, a constraint counting argument [6] shows that information derived from experiments with completely 
non-oriented particles is deficient for extracting 3D structure. In this paper we revive the prospects of Kam's 
idea in the 3D realm provided the particles in the ensemble can be partially oriented. The type of alignment 
required is illustrated with pyramidal particles in Figure 1. As a result of native or engineered interactions 
and geometry, we propose that the particles adsorb on a substrate in such a way that a particular oriented axis 
within the particle is aligned with the substrate normal with near 100% probabilityLj. Other means of achiev- 



1 With considerable additional effort the method we describe might be extended to allow for multiple modes of adsorption, in 
other words, when the dice are not as strongly loaded to land on one particular face as considered here. 




Figure 1: Partially oriented particle ensemble realized as dice with a strong tendency to land on 
one particular face. Alternatively, partial alignment might be realized in the manner of membrane 
proteins embedded in a lipid membrane. 



ing this degree of partial orientation would achieve the same goal. The concentration of their energy density 
makes substrates more efficient alignment mechanisms than laser fields, but this comes at the expense of 
introducing background and making certain ranges of the orientation angle inaccessible to the experiment. 
Membrane protein complexes, axially oriented in a synthetic lipid membrane, a liquid substrate, would be 
natural targets for this technique Q. 

The fluctuations that form the basis of Kam's idea are not due to photon number fluctuations but arise 
from different realizations of the particle orientations. In our case, these are limited to rotations about the 
substrate normal. If the particles undergo rotational diffusion, then the duration of each intensity measure- 
ment must be short compared to the diffusion time scale. When diffusion is absent, the sample substrate 
must be translated in order that a different set of particles, with different orientations, intercepts the beam. 
In either case, and from a signal perspective, it is always advantageous to focus the beam on the smallest 
number of particles N, if this can be done so without adversely affecting beam divergence and the integrated 
fluence. That is because, while the mean intensities recorded at two positions on the detector are indepen- 
dent of N, their covariance — the source of structure information — is diminished by 1/N. Even when the 
background signal at the two pixels is perfectly uncorrected, in being independent of N it will dominate the 
error in the covariance estimation for large N. There are other factors, however, that rule out small N even 
when this does not sacrifice the quality of the beam. Chief among these is particle damage, which must be 
kept at a low rate if subsequent measurements are taken on a single set of rotationally diffusing particles. 



In a nutshell, our reconstruction method involves three stages. In the data collection stage, intensity 
correlations are recorded for five continuous parameters. Four of these are the x and y positions of two 



pixels on the detector and the fifth is the tilt angle 9 of the substrate. This 5D data set is reduced to a 3D data 
set in the data reduction stage. The 3D data comprise complex numbers I m (r, z), or angular harmonics, 
giving the Fourier series expansion coefficients of the particle intensity on families of circles in a frame 
fixed to the particle. The data reduction also provides an internal consistency criterion, similar to principal 
component analysis, for establishing ranges for the radial (r) and longitudinal (z) frequency content, as well 
as the maximum harmonic order M that can be extracted from the data. The data reduction stage leaves 
M/2 phases in the intensity reconstruction undetermined. Finally, in the particle reconstruction stage, two 
phase problems are solved jointly that determine the 3D particle density from the (3D) angular harmonics. 
In addition to the usual phases encountered when reconstructing from intensity data, the M/2 unknown 
phases in the angular harmonics from the second stage are determined from the property that the intensity is 
consistent with a particle of compact support. 



2 Data collection 



The coordinate frame for data collection is defined by the incident x-ray beam, as z-axis, and orthogonal x 
and y axes in the plane of the detector. With the sample located at the origin, the unscattered beam would hit 
the detector at x = y = 0, z = Z. The axes on the detector are chosen with reference to the axis about which 
the substrate holding the sample is rotated. For tilt angle 9 = the substrate and detector planes are parallel 
and in-plane axes x' and y' on the substrate coincide with those on the detector. In all the measurements the 
substrate is tilted about the y' axis, leading to the following relationship between the lab and substrate axes: 

x = cos0x' + sin#z' (la) 

y = y' (lb) 

z = — sin#x' + cos#x'. (lc) 

The z' axis is the substrate normal and coincides with an alignment axis fixed in the body of each particle 
adsorbed on the substrate. A third set of axes, fixed in a particular particle, is free to rotate with respect to 
the substrate about its z' axis by some angle a: 
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x = cos ax + sin a y (2a) 

y = — sin a x" + cos a y" (2b) 

z' = z". (2c) 

In the experiment there is no control over the angle a and we assume it takes uniformly distributed random 
values for each particle. 

We can now express the photon wave vector change, seen in the laboratory, in the frame fixed to a 
particular particle. For intensity correlation measurements we are interested in the intensities at two points 
on the detector plane, that we label 1 and 2. The source of the scattered photons is the same particle, 
since the covariance vanishes for different particles whose rotation angles a are assumed to be independent. 
Using equations (Q} and (O, the wave vector transferred to photons detected at coordinates x\ and y\ on the 
detector is 

* - (£) * 

qi = (xi cos 9 cos a — y\ sin a) x" + {x\ cos 9 sin a + y\ cos a) y + (cci sin 9) z" (4) 



where A is the wavelength and there is a similar equation for q 2 . In subsequent equations we will express 
all wave vectors, as here, in units of 2ir/(\Z). In (@]) we made the small angle approximation, neglecting 
corrections of order (x/Z) 2 and (y/Z) 2 . 

For both pixels in our correlation measurement we define coordinates r, (j) and z as follows: 

x\ cos 9 = n sin </>i (5a) 

j/l = n cos 4>x (5b) 

xi sin# = zi, (5c) 

and a similar set of definitions for pixel 2. Let I(x, y; 9) be the (time integrated) intensity recorded in 
one measurement at detector pixel coordinates x and y and tilt angle 9. In the experiment we measure the 
following correlation function: 

C(x 1 ,y 1 ;x 2 ,y2;9) = cov (I(x 1 ,y 1 ;9),I(x 2 ,y 2 ;9)) . (6) 

The covariance is estimated in the usual way, from averages of intensities and their products taken from a 
set of diffraction patterns at fixed tilt. Let I p (x, y, z) be the intensity (squared Fourier magnitude) of the 
particle density as a function of the scaled wave vector components in the body-fixed coordinate system. 
Using dU) and the definitions ((5]>, the covariance is expressed in terms of the particle intensity as 

cov (I(x 1 ,y 1 ;9),I(x 2 ,y 2 ;9)) /N(9) = (7) 

( ip (n sin (</>i - a), n cos (0i -a),Zi) I p (r 2 sin (0 2 - a) , r 2 cos (cp 2 -a),z 2 )) a 

- (J p (nsin(0i -a),ricos(0i - a),Z\) } a (J p (r 2 sin(0 2 - a) , r 2 cos (<p 2 -a),z 2 )) a , 

where ( • • • } a is an average over the uniformly distributed particle rotation angle a and N(9) is the number 
of particles intercepted by the incident beam. In this equation the definition of the particle intensity I p has 
absorbed several factors, including beam fluence, that we assume are kept constant. In the following we 
further assume that the particle number density on the substrate is constant so that (again omitting a constant 
factor) we have 

N(9) = l/cos9. (8) 

Rotating both the tilted substrate and the detector by it about the beam axis should not change the corre- 
lations. This implies the following symmetry property of the correlation function: 

C(x 1 ,y 1 ;x 2 ,y 2 ; 9) = C(-xi, -yy, -x 2 , -y 2 ; -9). (9) 

It is therefore sufficient to only collect data for positive tilt angles. 



3 Data reduction 



As a first step in reducing the data we represent the particle density in terms of its angular harmonics with 
respect to rotations about the body axis parallel to the substrate normal: 

/ p (r sin </>,r cos </>,*) = ^e im ^/ m (r,z). (10) 



Using this definition in f7]) and the reality property of I p we can evaluate the average over a and obtain 

C(x 1 ,y 1 ;x 2 ,y2;9)cos0 = Y, eimi<t ' 1 ~ <t ' 2)l rn(ri,zi)I^(r2,Z2). (11) 

We can decouple these equations by averaging both sides of (fTTT) against the factor e~ im (^ 1_< ™ in such a 
way that 0i — 02 is uniformly sampled on an interval of 2tt radians. Moreover, the decoupling is especially 
simple if we can do this while keeping the parameters r\,Z\, r2 and Z2 all constant. 

The desired decoupling is achieved by a one-parameter average for which we use either 0i or $2, de- 
pending on the value of 

A=^. (12) 

T2Z1 

When \A\ < 1 we use 0i as the parameter; otherwise, by switching the labels 1 and 2 we restore \A\ < 1 
and the parameter becomes 02. In the following we assume the labels are such that \A\ < 1, making 0i the 
averaging parameter. Taking the ratio of equations (l5al) and (l5cl ) for both labels we obtain an equation for 
the second azimuthal angle in terms of the first, 

02(0i) = arcsin (^4sin0i), (13) 

where we take the branch that is continuous with 02 — » in the limit A — > 0. In the opposite limit, A — > ±1, 
02 is a sawtooth of amplitude tt/2. There is a second branch ir — 02 (0i) that we discuss below. For all A 
(in the range \A\ < 1), the function 

h) = 01 " 02(01 ) (14) 



is monotonic and ranges over exactly 27r radians whenever 0i does. Since we perform the average with 
respect to the parameter 0i, we need to include the Jacobian 



in order that </?(0i) is uniformly sampled. 
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t - AC ° S ^ (15) 

\/l - A 2 sin 2 0i 



Equations ([5]> and the companion equations involving label 2 determine how the detector parameters x\, 
Vi, %2, V2 and 6 need to be varied along with 0i and 02(0i) so that the body-fixed parameters r\, z\, T2 and 
Z2 are kept constant. Comparing d5al ) and ( f5cT > (and their companions) we see that the tilt angle must satisfy 

arctan I ) = Q{&\) = arctan I — — - ) , (16) 

Vrism0i/ V r 2sm0 2 (0i)y 

where compatibility follows from our definition ( fT3l of 02(0i). We take the branch of {<$>{) that lies in 
the range \9\ < tt/2, the center of which (6 = 0) corresponds to normal beam incidence on the substrate. 
Formulas for the other detector parameters are summarized below: 

(17a) 

(17b) 



Xi(0i) = 
£2(0l) = 

yiOi) = 
2/2(01) = 


sin#(0i) 

Z2 


sin#(0i) 

= 7"1 COS 01 

= r 2 cos0 2 (0i). 



(17c) 
(17d) 



Performing the one-parameter average of the correlation function we arrive at the decoupled equations: 

C m (ri,z 1 ;r 2 ,z 2 )= [ * ^ ./(^e"™^) coBO(<h.)C(xi(<h),Vl(<l>i)\ x 2 (0i),y 2 (0i);^(0i)) (18a) 
Jo 2tt 

= I m (n,z 1 )I* m (r2,z 2 )- (18b) 

Before we study equation (TT8l in greater depth we examine the second branch of the solution to (PT3l given 
by 

<P' 2 (<P 1 )=7T-MM- (19) 

Switching to the new parameter 

01 = 01 + 7T, (20) 

we find that the previously defined functions are changed in simple ways. Starting with 

2 (0i) = -0 2 (0i) (21) 

we find 

ip'ifa) = 4>i~ 2 (0l) = 0i " <M0i) - 2vr (22) 

= ^(0 , 1 )-2tt (23) 

so that (up to an irrelevant 27r shift) the averaging phase for the second branch is the same as for the first 
branch when expressed in terms of the parameter 4>'\ ■ The same is true of the Jacobian: 

J'{<h) = J(0i). (24) 

Finally, using (fT6l) and (fT71) one easily verifies the tilt and detector parameters for the second branch simply 
undergo a reversal of sign: 

(25a) 
(25b) 
(25c) 
(25d) 
(25e) 

The correlation average for the second branch is thus identical to (I18ab with all the arguments of the corre- 
lation function replaced by their negatives. The second branch therefore provides no additional information 
once symmetry (O is taken into account. 

In equation dl8a| ) we average the correlations C, between various pairs of pixels on the detector and 
various tilts, for particular harmonics m and define a correlation C m intrinsic to the particle in that its 
arguments are spatial frequencies r and z in the body-fixed frame. The two pixels move on circular arcs in 
the average, since equations © and their companions for pixel 2 imply 

xl + yl = r\ + z\ (26a) 

xl+yl = rl + zl, (26b) 

where the right sides are constant in any average. The left sides in these equations represent the magnitudes 
of the wave vector changes seen in the laboratory while the right sides are the magnitudes of the spatial 
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Figure 2: Pairs of detector pixels (labeled 0-12) in one averaging orbit (fT8l > for the case n = 6, 
•2i = 5, r2 = 9 and Z2 = —6. Lines through the circles show the tilt angle 0, where we have used 
symmetry © to keep the range of tilts between and tt/2. The tilt is smallest at pairs 3 and 9 and 
reaches 7t/2 (grazing incidence) at pairs 0, 6, 12. 



frequencies in the particle responsible for those changes. Figure 2 shows the orbits of two pixels and the 
corresponding tilt in an average for one particular choice of r±, Z\, T2 and z-i- Since all orbits visit tilt angle 
9 = it/2, where the beam is parallel to the substrate, the correlation function needs to be interpolated over 
the inaccessible range before the average (TT8T ) is evaluated. By using the symmetry (© it is only necessary 
to acquire data for 6 between and tt/2. 

In the final step of data reduction the angular harmonic amplitudes I m (r, z) are extracted from the aver- 
aged correlation function C m . For fixed m, the function C m (ri, Z\\r2-,z-i) is a Hermitian matrix with row 
indices n, Z\ and column indices r-i, z<i- The content of (118bl) is the fact that this matrix has rank one and 
its unique eigenvector with non-zero eigenvalue, up to an overall phase, is I m (r, z). 



The procedure for extracting I m (r, z) is similar to principal component analysis and tests the integrity 
of the data prior to the phase reconstructions that follow. First a set of radial (r) and longitudinal (z) spatial 
frequency samples are selected that, by d26l >. lie within the beamstop and edge of the detector. The density 



of the samples is set by the speckle scale of the single particle diffraction pattern. Next, the matrices C m 
are calculated by averaging correlations for the chosen r, z samples according to (1 18al) and for |m| < M/2, 
where the maximum harmonic M is set by the angular speckle scale at the edge of the detector. For each 
matrix C m one then obtains the dominant normalized eigenvector V m (r, z) and corresponding eigenvalue 
X m and forms the estimates 

I m (r,z) = \/X^V m (r,z). (27) 

From the relative magnitudes of the subdominant eigenvalues of C m one can assess the quality of the data at 
the chosen resolution. If these are significant, then it may be necessary to work with lower resolution blocks 
of C m , reduce M, or fix the interpolation of the correlation function at 9 = tt/2 until single-eigenvector 
dominance is realized for the C m . On the other hand, if no amount of data truncation satisfies this condition 
then the data are too flawed, e.g. by background, to proceed any further. 

We use the ratio of the largest eigenvalue magnitude to the sum of the magnitudes of all the eigenvalues 
of C m as an internal measure of the consistency of the reduced data, 

WC II 

ll^mlloo /oo\ 

which, as expressed by the notation, is the ratio of the spectral and trace norms of the matrix C m . In the 
absence of noise and other complications a m = 1. 

Since the phases of the eigenvectors V m in (127T ) are undetermined, our information about the intensity is 
incomplete at the level of M/2 unknown phase angles. These angles have to be reconstructed on the basis of 
additional constraints before the 3D intensity of the particle can be synthesized from its angular harmonics 
I m {r, z). Since M scales only with the diameter of the particle, and not its volume, this missing information 
is relatively minor in comparison with the many phases that must be reconstructed to obtain the particle 
density from the intensity. An algorithm for reconstructing both types of phases is described in section @] 

In contrast to the m ^ harmonics I m (r, z) that require correlation analysis, the m = harmonic is 
obtained by a straightforward average of the intensity and is therefore much easier to acquire. There is also 
no phase ambiguity in Io(r, z). Instead of © we define 

A{x,y;9) = & ve{I{x,y;9)) (29) 

and © is replaced by 

ave (I(x,y;9)) /N(9) = (J p (r sin(</> — a),r cos {<f> — a),z) ) a (30a) 

= I (r,z). (30b) 

We can aggregate the three parameter data on the left in a way that tests the two parameter model on the 
right. One method is to use the earlier parameterization (fTBT ) and (TTTT ) to define harmonics 

A m (r, z)= ^ ^ e - im * cos 9(d>) A (x (<£), y (</>); 9^)) (31) 

Jo 2?r 

such that Ao(r, z) becomes the estimate of Io(r,z) and power in the nonzero harmonics indicates data 
inconsistency. 



3.1 Data reduction for axial projections 

The problem of reconstructing an axial projection from data at tilt angle 8 = has been studied extensively 
by Dilano and co-workers [3, 4]. In this section we consider this special case and recover the main result of 
the previous studies. The final step of our data reduction, however, is different and is the starting point for 
the improved phasing method described in section 14.11 

By equations (f5]) the body-frame wave vector z vanishes at tilt angle 8 = 0. The corresponding limit of 
the particle intensity harmonic function I m (r, 0) provides information about the projected particle. Equation 
(fTTb reduces to 

C(x 1 ,y 1 ;x 2 ,y 2 ;0) = £ e im ^-^I m (ri,0)I^(r 2 ,0). (32) 



where the pixel coordinates are now simply 



^i(0i) = risin^i (33a) 

yi(4>i) = n cos 0i, (33b) 



and similarly for pixel 2. Specifying the angle of pixel 2 as </> 2 = <fii — <p, the right side of (1321 is independent 
of 4>i and we may average the data on the left side with respect to this angle. Decoupling with respect to the 
harmonic m is now accomplished by a simple integral with respect to ip: 

C m (ri,0;r 2 ,0) = y'^^e- im ^ /' 7r ^C'(s 1 (0i),yi(0i);x 2 (0i-^),y2(0i-^);O) (34a) 

= i m ( n , o)/;;(r 2 ,o). (34b) 

As in the general case, d34bb implies that the intensity harmonics are given by 

l m (r,0) = y/KMnir), (35) 

where V m (r) is, up to a phase, the unique eigenvector of the Hermitian matrix C m (n, ; r 2 , 0) with non- 
zero eigenvalue A m . The integrity of the data is assessed using the ratio a m (T28T ) as in the general case. 
The harmonic Io(r, 0) is just the z = 8 = limit of (I3TT) and corresponds, for m = 0, to a simple angular 
(powder) average of the intensity. 

We illustrate data reduction for the axial projection of a simple test particle whose intensity is shown in 
Figure 3. Uncorrelated, uniform amplitude and normally distributed noise r\ was added to the true intensity 
correlations C to simulate the effects of background. Our signal-to-noise ratio in this model is defined by 

SN = C rms /r/ rms , (36) 

where the root-mean-square amplitude of C is evaluated as a uniform average over all the pairs of measured 
pixels. 

After performing the angular averages (I34ab we extract the angular harmonics I m (r, 0) as the dominant 
eigenvectors of the matrices C m . The result of this for the test particle data with zero noise is shown 
in Figure 4. That the harmonics I m (r, 0) have a lower triangular support is a direct consequence of the 
uniform size of speckles in the intensity I p . Because our algorithm for extracting normalized eigenvectors 




Figure 3: Single-particle intensity I p (x, y, 0) of the test particle used in the reconstruction demon- 
strations. The dark circle at the center corresponds to the missing beam stop data. 

used an unspecified convention for the phase, there is no significance to the global phase of any of the rows 
of the harmonics plotted in Figure 4 (an arbitrary color-wheel rotation may be applied to each one). The true 
phases a m have to be determined in the final particle reconstruction stage. 

When noise is added to the intensity correlations the harmonics I m (r, 0) shown in Figure 4 become 
noisy as well. This is well diagnosed quantitatively through the principal eigenvector dominance measure 
a m and is shown plotted for our test particle in Figure 5 for three noise amplitudes. We will see that for our 
particular test particle as few as 10 reliably extracted harmonics (SN = 200) is nearly sufficient for density 
reconstruction while less than 5 harmonics (SN = 5) is insufficient. 



4 Particle reconstruction 



As remarked in connection with equation (|27T ). additional constraints are required to determine the M/2 
phases of the intensity angular harmonics I m (r, z). Non-negativity of the intensity is clearly one constraint 
that can be used. This constraint, however, is weak in comparison with the much stronger constraint that 
the intensity is the squared Fourier magnitude of a particle density having compact support. And since the 
particle density is our actual goal, we should approach the reconstruction of the M/2 intensity phases as 
part of a single process applied to intensity and density jointly. As we will show, the joint reconstruction 
is a relatively straightforward application of a general procedure once the constraints and the projections to 
them are set down (7J. 
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Figure 4: Result of the data reduction stage in the reconstruction of the axial projection of a test 
particle density. The plot shows the complex angular harmonics I m (r, 0) derived from the principal 
eigenvectors of the angularly averaged correlation matrix (I34ab . for each m, of the intensity in 
Figure 3. Magnitude and phase are rendered as brightness and hue, respectively. There is a global 
phase ambiguity for each m (row), the angles a m , that the data reduction does not determine. 

There are three constraints that apply to the pair {p, I}, where p is the Fourier transform of the particle 
density p and we have dropped the subscript p on the particle intensity /: 



support : 

data : 

compatibility : 



supp(p) C S, 

Jo* ^e~ tmip I(r Mi) v : r (os p. ; ) 



l I m {r,z), 



I. 



(37a) 
(37b) 
(37c) 



The first two apply to p and / individually and can therefore be combined into a single constraint projection. 
Projecting p to the support constraint is no different from how this is done in the standard phase reconstruc- 
tion problem. First p is Fourier transformed, the resulting density p is set to zero outside the support region 
S and also in its interior when the density is negative, and the result of this is inverse Fourier transformed to 
give the projection p' . The data constraint involves the harmonics I m (r, z) given by the data reduction and 
the unknown phases a m . Details of the corresponding constraint projection, which in effect determines the 
phases a m , are given below. 

With the support and data constraints combined, the compatibility constraint becomes the second con- 
straint in the general scheme where a problem is divided into just two easy constraints. This division of 
the reconstruction problem is similar to the strategy known as divide and concur |8). The variables p and 
I are to a large extent redundant representations of the same information and by making them independent 
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Figure 5: Dominance of the principal eigenvector as measured by the quantity a r , 
from the test particle data reduction and three level of noise. 



defined in (l28l) 



we facilitate the projection to the "divided" constraints (support and data). In our case "concurrence" of the 
redundant variables is imposed by the compatibility constraint, which takes a nonlinear form. 

With any projection to a particular constraint set comes the understanding of the distance in the space 
of variables that the projection minimizes. In our case the variables are of two types: complex Fourier 
amplitudes p(x, y, z) and non-negative intensities I(x, y, z), both sampled on the same Cartesian grid. The 
Euclidean distance in this space of variables, between reconstructions {p, 1} and {/>', I'}, should be written 

^-^ w *-^ 

x,y,z x,y,z 

where the parameter w addresses the fact that p and / first of all have different units. Since the compatibility 
constraint satisfied by the projected variables, at each grid point (x, y, z), is the equation 

\p{x,y,z)\ 2 = I'(x,y,z), (39) 

we see that w has the same units as /. Although the precise value of w has a direct effect on the projection 
{p, 1} — > {/3', I'}, the projected values will always satisfy d39l ), which is independent of w. The performance 
of the reconstruction algorithm, on the other hand, can be adversely affected when w is either very small 
or very large. In those limits one of the variable types, p or /, will be much more compliant than the 
other and as a result the constraint associated with it (support or data) will be explored more rapidly than 
the constraint imposed on its stiffer counterpart. We obtained good results in the numerical experiments 
described in section |4~T1 when w was set to the maximum measured value of Io(r, z), near the minimum of 
q 2 = r 2 + z 2 , and making w larger or smaller than this by a factor of 2 had little effect on convergence. 
If necessary w can be allowed to vary with q while only slightly complicating the projection to the data 
constraint, the only other projection dependent on w. 



Projecting to the compatibility constraint ( [391 reduces to a one-parameter numerical optimization at each 



grid point (x, y, z). Given input variables p(x, y, z) 



ve 



and I(x, y, z) = I, the output of the projection 
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is given by p'(x, y, z) = v'e lt and I'(x, y, z) = I 1 where v' and I' satisfy 

v' 2 = I' (40) 

and minimize 

(v' - vf + (/' - if/w. (41) 

The resulting cubic equation for v' always has a unique, positive solution which can be tabulated in advance 
for efficiency 

The nontrivial part of projecting to the data constraint involves only those spatial frequency combinations 
(r, z) whose magnitude q lies between q min at the beamstop and q max at the edge of the detector. For 
1 < <Zmin, the projection merely copies the input intensities into the output intensities, while for q > q max 
(in the corners of the intensity grid) the intensities are set to zero. 

We uniformly sample the half-annulus defined by q m - m < q < q ma , x and r > on a Cartesian grid with 
equal spacings for r and z, and finer than the speckle scale. Every sample (r^, Z{) defines a circle of radius 
ri in a plane at level Zi of the intensity. The first step of the (non-trivial) data projection is extracting the 
angular harmonics of the input intensity by Fourier transforming the intensity on these circles: 



M-l 

Imi = j-7 ^2 e" Jm ^ I(n sin (fj , ^ cos <pj ,Zi), (42) 

where the angular samples (fj = 2-kj/M are equally spaced and M is chosen to match the angular speckle 
scale at r = q max . From the data reduction stage we know these harmonics up to an overall phase: 

All that remains is to determine the angles a m by minimizing 

M-l A/-1 

/ _, / _, \ 1 mi 1 mi\ — / J / ^ | e 1 mt 1 mi\ ■ V^+J 

m=0 i m=0 i 



The metric of (|44l) is equivalent to the metric (1381 ) by Parseval's theorem. Since ao = 0, the projection 
simply replaces Iq{ by the powder average 1^. For m/0 the minimization of (l44l yields 



g\S2I; m i m X (45b) 



oc m = ar, 



So far we have defined two projections 

P D ,P C : {p,I}^{p',I'}, (46) 

where Pp, the "divided constraint projection", imposes the support constraint (I37ab on p and the data con- 
straint (I37bb on /, while the "concurrence projection" Pq imposes the compatibility constraint (I37cb between 
p and / (while ignoring the support and data constraints). As distance minimizing operations, these are in 
a sense "greedy" moves to their corresponding constraint sets. Solutions {/> S oi,^soi} to the reconstruction 
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Figure 6: Iterations 40, 70 and 200 in the simultaneous reconstruction of intensity (top) and density 
(bottom) of the test particle axial projection (the greek letter a) using the algorithm described in 
section 14.11 The square frames in the bottom row are the boundary of the support. 



problem have the property that they are fixed by both projections, and in particular, p so i is simultaneously 
consistent with the intensity correlation data and a compact support. There are provably convergent algo- 
rithms for finding fixed points in a general projection setting when the constraint sets are convex. Because in 
our case constraints (137bl ) and (137cb are non-convex, convergence cannot be proven and we are forced to use 
algorithms that have a demonstrated record of discovering solutions even for difficult non-convex problems. 
We use the difference map iteration Q 



{p, 1} -> {p, 1} + P C (2P D {p, 1} - {p, I}) - Pd{p, I}, 



(47) 



which is equivalent to Fienup's hybrid input-output iteration [9] in the standard reconstruction problem with 
support and Fourier magnitude as the two constraints (Pc — > Ps, Pd — > Pf)- For our initial iterates we 
generate a random positive density p (uniformly distributed contrast) on the support and from this obtain p 
and I = \p\ 2 , a configuration that only satisfies the compatibility constraint. Iterations are performed until 



the metric A of the updates (1381 ) fluctuates about a steady state and shows no sign of further decrease. 



4.1 Reconstruction of axial projections 

The specialization of the reconstruction algorithm to axial projections amounts to simply setting z = in 
the equations of the previous section. In the flat Ewald sphere limit we can also take advantage of Friedel 
symmetry, thereby reducing by a factor of 2 the work in projecting to the data constraint (the angular Fourier 
transforms act on a periodic function on the domain < ip < ir). 

The advantage of simultaneously reconstructing intensity and particle density is evident from the frame 
by frame development of these in the early iterations of the difference map. There is a "locking" of the 
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Figure 7: Time series of the difference map error metric A in test particle axial projection recon- 
structions for three values of noise. The corresponding particle densities are shown in Figure 8. 

angles a m , that determine the intensity reconstruction, at about the same point in time that the density is 
well localized in the support, a necessary condition to have speckles in the intensity. In the subsequent 
refinement iterations the particle density and intensity often executed significant (synchronized) rotation. 
We deliberately chose a rather loose support to enable this freedom in the reconstruction. Figure 6 shows 
three frames in a reconstruction with low noise (SN = 10 4 ). The difference map error metric is plotted 
in Figure 7 for three levels of noise and not surprisingly shows that the fixed-point property (A — > 0) is 
strongly compromised when the noise is high. The corresponding reconstructions are shown in Figure 8. 

Given the relatively small number of phase angles required for intensity reconstruction, the weaker con- 
straint arising from the support of the intensity Fourier transform — the density autocorrelation — might 
avoid the complications in the simultaneous intensity/density reconstruction described above. The latter 
constraint is a weaker constraint on the intensity in that the Fourier transform is only required to have a par- 
ticular support, and not necessarily the structure of an actual autocorrelation. We tested this simplification 






Figure 8: Reconstructions of the test particle at SN values 5, 200, and 10 4 . 
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Figure 9: Difference map error metric when reconstructing just the test particle intensity from low- 
noise data and the autocorrelation support constraint. Convergence is much slower with this weaker 
prior constraint than for the joint intensity/density reconstructions (Fig. 7). A well reconstructed 
intensity (not shown) first appears only after 500 iterations. 



for our test particle by pairing the projection to the data constraint (I37bl ) with a simple autocorrelation sup- 
port projection in the difference map scheme. This completely avoids manipulating the particle density and 
the non-linear compatibility constraint (I37cb . Whereas this simpler approach did succeed in reconstructing 
the intensity, the progress of the algorithm, shown in Figure 9, was much slower and more sensitive to noise. 
And since the density is always the final goal of the reconstruction anyway, this two-stage approach does 
not appear to offer any advantages over the simultaneous intensity/density reconstruction. 
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